Terahertz relativistic spatial solitons in doped graphene metamaterials 
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We propose an electrically tunable graphene-based metamaterial showing a large nonlinear optical 
response at THz frequencies, which we calculate analytically for the first time to our knowledge and 
arises from the intraband current. The structure sustains a novel type of stable two-dimensional 
spatial solitary wave, a relativistic version of the Townes soliton. These results can be also applied 
to any material exhibiting a conical dispersion with massless Dirac fermions. 
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Introduction — Graphene is a two-dimensional, one- 
atom thick allotrope of carbon that has been in the spot- 
light since its experimental discovery and isolation in 
2004 [1| , and can be considered a unifying bridge between 
low-energy condensed matter physics and quantum field 
theory, as its two-dimensional quasi-electrons behave like 
massless "relativistic" Dirac fermions, very similarly to 
electrically charged 'neutrinos' [2|, |3[. 

Graphene holds the promise for building advanced 
nano-electronic devices, due to its unconventional elec- 
tronic properties [3|. Furthermore, it also exhibits very 
unique optical properties, especially in the terahertz 
(THz) frequency range. To date, novel photonic devices, 
such as THz devices [^ , optical modulators [5] , photode- 
tectors [65 and polarizers [7^ were successfully realized. 

In recent years, the huge and largely unexplored po- 
tential of graphene for nonlinear optical applications has 
been outlined. An extremely strong nonlinear optical re- 
sponse in the THz regime has been investigated [Sl, l9|. 
Preliminary experimental results include ultrafast sat- 
urable absorption 1 101 . Ill] and the observation of strong 
four- wave mixing |12|, which are the building blocks of 
nonlinear optics [13]. Specifically, four- wave mixing im- 
plies the existence of modulational instability and optical 
solitons in graphene, a significant topic that has not been 
previously investigated. 

In this paper, we follow the footsteps of a series of 
seminal papers by Mikhailov (see Refs. [14]) based on 
the semiclassical kinetic theory, and we derive analyti- 
cally the intraband optical current of a doped layer of 
graphene. We prove that this theory is also consistent 
with the more precise quantum approach of Ishikawa [9] , 
based on the Bloch equations derived from the single- 
electron Dirac equation. For excitation frequencies in 
the THz gap, and neglecting the interband transitions 
(valid for photons below the Fermi energy of doped lay- 
ers), the two approaches give the same result for the in- 
traband current. We apply the above results to describe 
self-focusing of two dimensional Townes-like solitons in 
an electrically tunable metamaterial made of several lay- 
ers of doped graphene, interspaced by layers of silica and 
silicon with thickness much smaller than the wavelength 
of the incident light. 




FIG. 1: (a) Graphene conical dispersion with doping. Intra- 
band and interband optical transitions are shown, (b) Geom- 
etry of the proposed multilayer metamaterial. The structure 
is made of graphene-silica-silicon layers, with total thickness 
L, much smaller than the wavelength of the THz beam. Each 
layer of graphene is doped by a gate voltage Vg. 



Background — We consider an electrically doped 
graphene system with a positive gate voltage Vg. As 
shown in Fig. [TJi, the electron energy dispersion in 
the conduction band is given by the Dirac spectrum, 

ep = ^fP, where p = |p| = 



+ Py is the total momen- 
tum, p = {pxiPy)i and '^F — c/300 if the Fermi velocity 
with c the vacuum light speed. The Fermi energy eF can 
be largely controlled by the voltage Vg perpendicularly 
applied to the graphene-Si02-Si multilayer (Fig. [TJd). 
The velocity operator for the quasi-electrons is given by 
V = Vpep = vfp/p- The electron momentum distribu- 
tion /p(r,t) in the collisionless approximation is solution 
of the Boltzmann-Vlasov kinetic equation: 

ajp(r, i) + V V/p(r, t)+¥- Vp/p(r, t) = 0, (1) 

where r and t are respectively space and time co- 
ordinates, and F = — eE is the force due to the 
electric field E(r,t) (e > here and in the follow- 
ing). Due to Jeans' theorem [15|, any function of the 
constants of motion is a solution of the Boltzmann- 
Vlasov equation, and assuming for simplicity homo- 
geneity in the {x^y) plane (i.e. V/p = 0), an ex- 
act solution of Eq. ([!]) is the Fermi-Dirac distribu- 
tion at temperature T for negligible interband tran- 



sitions, namely /p(t) = Tt [Px - Po,x{t),Py - Po,y{t)], 
where Tt [Px.Py] = [1 + exp {(ep - ep) /(/cbT)}]"\ ks is 
the Boltzmann constant, po(^) = —eA{t) is the electron 
momentum transferred by the radiation field, A(t) = 
— J'E{t)dt is the vector potential, and ep corresponds 
to an electron surface density n^ = e^/{h'^v^7r). 

Calculation of total intraband current — The total 
electric current is given by J = —rf^fe-e / v/pdp = 

~"(I^^'^f/(p/p)^t[p - Po(^)]c^P, where ^s = 2 and 
^v = 2 are respectively the spin and valley degeneracy 
factors. This gives the intraband current that is re- 
sponsible for the strong THz nonlinear ity of graphene, 
as demonstrated by Mikhailov 12|, [l^ . 

More recently, Ishikawa [9] introduced Bloch-like equa- 
tions deduced from the one-electron Dirac equation. In 
his formalism, he starts from the Weyl equation for the 
charged neutrino, icr^d^il) = 0, where d^ = (%^9t, V^) 
is the pseudorelativistic derivative, Vx = (dx^dy)^ and 
a^ = (o-^^cr) is the Pauli matrices vector, and cr^ is 
the 2x2 identity matrix. Expanding into space and 
time variables, this translates into the wave equation 
for the electronic spinor with momentum p, which reads 
ihdtil^p = 'UF(cr • p)'0 = Hoipp^ with the momentum op- 
erator p = —ih\/±. The interacting theory is directly 
implemented via the minimal substitution: ihdtil^p = 
vy {(J - [p + ^A(t)]) V^p. Without loss of generality, we 
can assume that the radiation field is polarized linearly 



along the x direction, parallel to the graphene plane. The 
Weyl equation for one electron is therefore 



ihdt^l^p — '^F 
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(2) 

From Eq. ([2]), one can construct the single-electron cur- 
rent Ji,p = —ev-pipj^aipp for a given electronic momentum 
p. Ishikawa [9] demonstrated that, for intraband transi- 
tions only (n = — 1 and p = in the notation of Ref. [9[), 
the single-electron intraband current reduce to 



Ji,p 



-evF 



{px^eAf^pl 



-.{Px^eA.py), 



(3) 



and J = (-|^^')2 j *^i,p^T[Px-,Py]dY>' This result agrees 
with the Boltzmann- Vlasov equation approach given 
above. It is important to note that, following Refs. |14| . 
the intraband current dominates the interband current 
for photon energies hw < ep and for /cbT <C ep. 

In this paper, we report the first analytically calcu- 
lation of the total macroscopic intraband current of a 
single graphene layer of thickness d at low temperature 
(T ^ 0). The final result for the x component (i.e. the 
only non-vanishing component of the two-dimensional 
current) is given by 
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where we have defined the elliptic integrals E±{x) = 
£±{^\x), with £±iO\x) = J^{l-xsm^O)^^/^dO. 

By defining the dimensionless variable ip = eA/p^^ 
where pf = ep/'^F is the Fermi momentum, and elec- 
tric monochromatic fields are scaled with the reference 



field ^0 = ujeF/{vFe). Fig. [2fa) shows the dimension- 
less quantity J2d('^) = -^2d/Jf (blue solid line) [where 
JF = —cvfPfK^^^) is the elementary Fermi current] 
which in terms of the ip variable is given by 
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showing the strong intrinsic nonlinear behavior of 
graphene layers when excited by THz radiation. Inter- 
estingly, J2d('^) can be roughly approximated in sev- 
eral ways, depending on the specific application. For 
instance, a hyperbolic tangent function approximation, 
J2d('^) — tanh('0), is useful for several estimates, for 



example in the calculation of the linear intraband con- 
ductivity, giving the correct aintra = e^ cf / {"Kh^ uj) with 
uj being the radiation frequency. We note that the lat- 
ter function has exactly the same asymptotic behavior of 
J2D for '^ ^ and ijj ^ oo. However the approximation 
might fail for more precise estimates around the most 
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FIG. 2: (a) Plot of the current J2d('0) (blue solid line), its 
hyperbolic tangent approximation (black dashed-dotted line) 
and its relativistic approximation (red dashed line). (b,c,d) 
Plots of analytical intraband current J2V>{t) when '0(t) = 
V^osech {t/to) cos(5t/to) (blue solid line), for '0o = 0.2, 1 and 3 
respectively, corresponding to the three black dots in (a). The 
lines for the tanh (black dashed-dotted) and the relativistic 
(red dashed) approximations are also shown. 



nonlinear region, namely ip = 1. Here, however, we pre- 
fer to use the less precise but more tractable approxima- 
tion J2d('^) — '^/yl + V^? which shows also the pseudo- 
relativist ic nature of the optical nonlinear ity treated here. 
This expression will allow us to treat the transition from 
real fields to envelopes in a straightforward way in the 
framework of the paraxial approximation, since the Tay- 
lor expansion of Eq. (|5j) does not work well due to its 
fictitious singularity at t/^ = 0. 

In Fig. [2l^a) J2d('0) and its two approximated versions 
are shown for comparison. In Figs. [2fb,c,d) we show the 
analytically calculated current J2d(^) for an example of 
pulsed excitation ^l)(t) = V^osech (t/to) cos(5t/to), where 
to is the pulse width, for different values of the light field 
amplitude '^o, showing the strong nonlinear temporal de- 
pendence of the intraband current. Curves obtained by 
using the tanh and the relativistic approximations above 
are also shown for comparison. 

A quick estimate of the third-order susceptibility Xgr is 
given by expanding J2d in powers of ip up to the third or- 
der: J2D — ip — ip^ /S, giving the nonlinear third order in- 
traband current J^^{E) = [eE/{2ujpF)f[evFpl/{7Th^)] = 
ujcoXgr E^' The order of magnitude of such susceptibil- 
ity for a monochromatic wave is thus given by Xgr = 
e^v^/iSneoh'^uj^eFd) - lO^^lO^'^xlnLa' depending on the 



intraband nonlinearity rapidly decreases when increasing 
the frequency. This is consistent with the estimates given 
by Mikhailov [12|, [l^. Such estimates, although relevant 
to retrieve the orders of magnitude involved, do not cap- 
ture the full complexity of the nonlinearity, which is not 
a simple Kerr nonlinearity. However, such a large third 
order coefficient places the graphene nonlinearity in the 
same category of the resonant nonlinear effects, such as 
two- level systems 16| or the excitonic nonlinearity [17|, 
but with the great advantage that the bandstructure of 
graphene is always resonant to optical excitations. 

Wave equation for a single graphene layer — The 
equation that regulates light propagation in presence 
of a single graphene layer is the conventional macro- 
scopic wave equation c^eoDA = J3d(^), where D = 
(es/c^)9f — V^, with eg is the substrate dielectric function 
at the selected frequency cj, and A(r,t) is the 3D vector 
potential. The current circulates in a very thin layer of 
thickness d c:^ 0.34 nm [18']. Thus we can model the 3D 
current with a rectangular function, with the single layer 
centered at z = 0: J3D(r, t) = J2B{x,y^t)R{z)/d, where 
R{z) = {sgn{z^d/2)^sgn{d/2- z)}/2, where sgn(x) 
is the sign function, normalized like J_ R{z)dz = d. 

Average medium theory of graphene metamaterial — 
In order to observe THz spatial solitons, we consider a 
doped graphene metamaterial as shown in Fig. 1(b). The 
system that we propose is a periodic multilayer based 
on graphene-silica-silicon layers, with total thickness L 
[see Fig. [1] (b)] where Si02 and Si are both transpar- 
ent for THz and optical frequencies [19]. Each graphene 
layer is doped with an electronic density n^ by applying 
a gate voltage Vg. The size L of the elementary cell is 
assumed to be much smaller than the incident monochro- 
matic wavelength, L <C A. This means that we can use 
an average medium approach by expanding the dimen- 
sionless vector potential in its Fourier components ^^n- 
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(6) 

where ko is the linear wave number, and z is the direction 
perpendicular to the layers [see Fig. [Hb)]. By retain- 
ing only the fundamental order in the Fourier expansion 
(with (j) = 00 ) 7 and after using the paraxial approxima- 
tion for reducing the order of the z derivative, one ob- 
tains: 
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iJ^^^,^R{z)dz = d/L. 



We now use the relativistic approximation of the full cur- 
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FIG. 3: (a) Soliton profiles of the fundamental mode (n — 0) 
and two higher- order modes (n = 1 and n = 2) as a function 
of dimensionless radius i?, for q — —0.7. (b) Same as (a), 
but for q — —0.1. (c,d) 3D plots of the fundamental and the 
2nd-order soliton of (a,b) respectively. 



graphene analytically, which dominates the electron dy- 
namics for THz excitations. Finally, stable Townes-like 
spatial solitary waves have been found to propagate in 
the longitudinal direction for realistic parameters. These 
results pave the way of a more extensive analysis of the 
mathematical structure and the physical content of the 
graphene Bloch equations for useful nonlinear optical ap- 
plications. Our theoretical approach is not restricted to 
graphene, but can be applied to all materials exhibiting 
a conical dispersion supporting massless Dirac fermions 
(for instance HgTe [26,]). 
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rent, namely J2d(0) — 0/[\/l + |^p/2], inside Eq. ([7]). 

Paraxial model and soliton solutions — After intro- 
ducing the scalings {x^y) = xo(X, y), z = zqZ^ xq = 
[7^fi^c^eoesL/(e^eF)]^/^ zo = 2/coxg, ko = ^Je^ujlc and 

T] 



^/>/2, one obtains the paraxial equation: 
idzT] + {d\ + dl)r] 



VT 



l^l' 



(8) 



For instance, for typical parameters cj/(27r) = 20 THz, 
Cs '^^ 4.5 [20], L :^ 2 /im and n^c^h^ 10^^ cm~^, one has 
ep — 259 meV, xq — 15 /im and zq — 415 /im. The 
scaling for the electric field is E^ = ujeYJiyYe) o:^ 30 
MW/cm^. For the above parameters the room temper- 
ature is a good approximation of the above results ob- 
tained for T = 0, since /cbT :^ 25 meV <C ep- 

Passing to cylindrical coordinates one must solve the 
following ODE: 



d^rjjR) 1 dr]{R) 
dR^ ^RdR 



1 



yrr^PF 



v{R) = 0, (9) 



where g' is a nonlinear wave number, and R = \JX^ + Y^ 
is the dimensionless radius, in units of xq. Solutions of 
Eq. (j9]) are Townes-like solitons (see Refs. [2l'-^) with a 
rather unconventional relativistic nonlinearity, which are 
stable in the sense of the Vakhitov-Kolokolov criterion 
due to the saturable type of the nonlinearity [25]. Some 
fundamental and higher-order soliton profiles are shown 
in Fig. [3] for different values of q. 

Conclusions — We proposed an electrically tunable 
metamaterial based on graphene-silicon-silica multilay- 
ers. We calculated the intraband current of doped 
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